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S ■ ABSTRACT 

^ , Cold accretion disks with temperatures below ~ 3000K are likely to be 

^ ', composed of highly neutral gas. The magnetorotational instability may cease 

CN ! to operate in such disks, so it is of interest to consider purely hydrodynamic 

^ i mechanisms of generating turbulence and angular momentum transport. With 

^ i this motivation, we investigate the growth of hydrodynamic perturbations in a 

Q\ ', linear shear flow sandwiched between two parallel walls. The unperturbed flow is 



(N 



similar to plane Couette flow but with a Coriolis force included. Although there 
are no exponentially growing eigenmodes in this system, nevertheless, because of 
^ I the non-normal nature of the eigenmodes, it is possible to have a large transient 

^ I growth in the energy of perturbations. For a constant angular momentum 

disk, we find that the perturbation with maximum growth is axisymmetric 
O I with vertical structure. The energy grows by more than a factor of 100 for 

I a Reynolds number R = 300 and more than a factor of 1000 for R = 1000. 

^ ' Turbulence can be easily excited in such a disk, as found in previous numerical 

simulations. For a Keplerian disk, on the other hand, similar perturbations 
^ I with vertical structure grow by no more than a factor of 4, explaining why the 

I same simulations did not find turbulence in this system. However, certain other 

two-dimensional perturbations with no vertical structure do exhibit modest 
growth. For the optimum two-dimensional perturbation, the energy grows by 
a factor of ~ 100 for R ~ 10^-^ and by a factor of 1000 for R ~ 10^. Such 
large Reynolds numbers are hard to achieve in numerical simulations and so the 
nonlinear development of these kinds of perturbations are only beginning to be 
investigated. It is conceivable that these nearly two-dimensional disturbances 
might lead to self-sustained three-dimensional turbulence, though this remains 
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to be demonstrated. The Reynolds numbers of cold astrophysical disks are 
much larger even than 10^; therefore, hydrodynamic turbulence may be possible 
in disks through transient growth. 

Subject headings: accretion, accretion disk — hydrodynamics — turbulence — 
instabilities 



1. Introduction 

The origin of hydrodynamic turbulence is not fully understood. Many efforts have been 
devoted to this problem, beginning with the early work of Kelvin, Rayleigh and Reynolds 
toward the end of the nineteenth century. However, despite a large number of investigations 
over the decades, the key physics is still poorly understood. One of the reasons for this 
is that there is a significant mismatch between the predictions of linear stability theory 
and experimental data. For example, plane Poiseuille flow is known to become turbulent 
in the laboratory at a Reynolds number R ~ 1000, whereas theory predicts that the flow 
is linearly stable up to i? = 5772. An even more severe discrepancy, one that is of direct 
interest to astrophysics, occurs in the case of plane Coucttc flow. Laboratory experiments 
and numerical simulations show that this flow can become turbulent for R as small as 
~ 350. However, theoretical analysis shows that the flow is linearly stable for all R up 
to infinity. Such a large discrepancy indicates that linear stability analysis, based on 
eigenspectra, is not the best tool for understanding the onset of turbulence. In this paper, 
we pursue a different approach, the so-called bypass mechanism to transition, which has 
been popular in the fiuid mechanics literature (e.g., Farrell 1988; Butler & Farrell 1992; 
Rcddy & Henningson 1993; Trefethen et al. 1993). We use this approach to study a possible 
route to turbulence in astrophysical disks. 

Accretion disks in astrophysics operate by transferring angular momentum outward 
by an effective "viscosity." Microscopic molecular viscosity is completely negligible, so it 
was recognized more than three decades ago (Shakura & Sunyaev 1973; Lynden-Bell & 
Pringlc 1974) that angular momentum transfer must occur via turbulence of some sort. 
However, the physical origin of the turbulence was not identified until the important 
work of Balbus & Hawley (1991) who identified the Magneto- Rotational-Instability (MRI) 
(originally discovered by Velikhov 1959; Chandrasekhar 1960) and showed that this linear 
instability will operate in the presence of very weak magnetic fields and will lead to 
magnetohydrodynamic (MHD) turbulence. The MRI is now accepted as the origin of 
turbulence in most accretion disks. Hawley, Gammie & Balbus (1995) showed that the MRI 



-3- 



dies out if the Lorentz force is turned off, and Hawley, Gammie & Balbus (1996) showed that 
the magnetic field dies oTit when the CorioHs force is turned off while retaining the Lorentz 
forces. In both of these situations, MHD turbulence is absent. In two subsequent papers, 
Balbus, Hawley & Stone (1996) and Hawley, Balbus & Winters (1999) showed through 
numerical simulations that, whereas pure hydrodynamic turbulence is easily triggered in 
plane Couette flow (which was already known) and in a constant angular momentum disk, 
turbulence does not develop in an unmagnetized Keplerian disk even in the presence of 
large initial perturbations. The authors argued on this basis that hydrodynamic turbulence 
cannot contribute to viscosity in accretion disks. 

Despite the above important work, there is reason to study hydrodynamic turbulence 
in astrophysical disks. Several accretion systems are known in which the gas is cold and 
largely neutral so that the gas and the magnetic field are poorly coupled. The MRI then 
becomes weak or may even cease to operate. In a series of experiments, Hawley, Gammie & 
Balbus (1996) and Fleming, Stone & Hawley (2000) showed that for a magnetic Reynolds 
number, Rm ~ 10^, the magnetic field is depressed and for Rm — 2x 10^ the magnetic field 
dies out. Thus, for Rm < 10^, MHD turbulence and the associated angular momentum 
transport switch off. Examples of systems in this regime include accretion disks around 
quiescent cataclysmic variables (Gammie & Menou 1998; Menou 2000), proto-planetary 
and star-forming disks (Blaes & Balbus 1994; Gammie 1996; Fromang, Terquem & Balbus 
2002), and the outer regions of disks in active galactic nuclei (Menou & Quataert 2001; 
Goodman 2003). Note that the magnetic Reynolds number may not be the only relevant 
parameter that determines the strength of the MRI; the magnetic Prandtl number may 
play a role (Gammie & Menou 1998), as well as ambipolar diffusion (Blaes & Balbus 
1994; Menou & Quataert 2001). Nevertheless, it seems reasonable to assume that, in cold 
astrophysical disks, the MRI will be sluggish or even absent. The question then arises: 
How can these systems sustain mass transfer in the absence of the MRI? What drives their 
turbulence? 

A number of ideas have been discussed in the literature in answer to this question. 
Gammie (1996) argued that the surface layers of cold protostellar disks would be ionized 
by cosmic rays and that this would enable accretion to proceed via the MRI within these 
layers. Menou (2000) suggested that angular momentum transport in quiescent cataclysmic 
variables could be induced by tidal perturbations from the binary companion star and 
showed that this mechanism could explain the occurrence of longer outburst recurrence 
times in systems with large binary mass ratios. In the case of the outer disks of active 
galactic miclei. Menou & Quataert (2001) showed that the MRI-stable regions are nearly 
always gravitationally unstable, so that the latter might drive angular momentum transport. 
A problem with these ideas is that each is invoked specifically for a particular class of 
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systems. While there is nothing in principle wrong with this, one wonders whether there 
may not be some more general mechanism for generating turbulent transport in MRI-stable 
astrophysical disks. 

Recent laboratory experiments on rotating Couette flow in the narrow gap limit with 
hnearly stable rotational angular velocity profiles (similar to Keplerian disks) seem to 
indicate that turbulence does manage to develop in such flows (Richard & Zahn 1999). 
Longaretti (2002) points out that the absence of turbulence in the simulations by Balbus, 
Hawlcy & Stone (1996) and Hawlcy, Balbus & Winters (1999) may be because of their small 
effective Reynolds number. Also, Bcch & Andcrsson (1997) have shown that turbulence 
persists in numerical simulations of sub-critical rotating flows, provided the Reynolds 
number is very high. Moreover, as already mentioned, it is well known that linear stability 
is no guarantee that a flow (whether rotating or not) will avoid becoming turbulent (for a 
detailed discussion see e.g. Swinney & GoUub 1981; Drazin & Reid 1983). In fact, since 
the celebrated work of Orr (1907), it has been known that linearly stable flows can exhibit 
significant transient growth in energy for certain initial perturbations. This fact provides a 
possible solution to the problem of explaining hydrodynamic turbulence in linearly stable 
systems. The idea is that the transient growth may allow perturbations to grow to a 
non-linear state, after which a sub-critical transition to turbulence may take place. This 
is called the bypass mechanism to turbulence. In the astrophysical literature, an early 
application of transient growth may be found in Goldreich & Lynden-Bell (1965; see also 
Goldreich & Tremaine 1978, 1979). 

The physics of transient growth has been discussed by a number of authors (Farrell 
1988; Butler & Farrell 1992; Reddy & Henningson 1993; Trefethen et al. 1993), who have 
shown that the growth results from the non-normal nature of the associated operator. The 
eigenfunctions of the linearly perturbed system are not orthogonal but are close to linearly 
dependent in nature, and as a result certain linear combinations of the eigenfunctions 
that are arranged to nearly cancel initially may develop considerable amplitude at later 
time when the degree of cancellation is reduced. Therefore, even in the absence of any 
exponentially growing eigenfunctions, the system is still able to exhibit transient growth. 
This idea has been discussed in the fluid mechanics literature for a number of years but 
has only recently been applied to astrophysical accretion disks. loannaou & Kakouris 
(2001) studied the global behavior of perturbations in an accretion disk, Chagelishvih et al. 
(2003) analysed a local 2-dimensional patch in a disk using a shearing box approximation 
and showed that strong growth is possible, and Tevzadze et al. (2003) showed that 
3-dimensional perturbations also undergo substantial transient growth, provided the vertical 
scale remains of the order of the azimuthal scale. Umurhan & Regev (2004) studied the 
non-linear development of the Chagelishvih et al. (2003) growing mode and Yecko (2004) 
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studied rotating shearing flows between walls. Recently, Johnson & Gammie (2005a,b) 
studied the evolution of a plane-wave type perturbation in thin low-ionization disks. 

The aim of the present study is to further explore the physics of transient non-normal 
growth of perturbations in cold accretion disks, with a view to understanding whether 
such growth could lead to hydrodynamic turbulence. Our aim is to present the analysis 
in such a manner that even readers from other branch of astrophysics (not familiar with 
the conventional fluid dynamical approach) will be able to follow and reproduce the results. 
Along with a companion paper (Afshordi, Mukhopadhyay & Narayan 2005; hereafter 
AMN05), we study both Keplerian and constant angular momentum disks, as well as 
plane Couette flow. Both papers concentrate on identifying the parameter regimes over 
which a large transient growth in energy is possible and studying the nature of the growing 
perturbations. While the present paper focuses on an eigenvalue analysis in Eulerian 
coordinates of flow between walls, AMN05 presents a Lagrangian analysis of an infinite 
shear flow. 

The plan of the paper is as follows. In §2, we present our basic model, beginning 
with a description of the equilibrium flow, then discussing the perturbation equations and 
eigenfunctions, and introducing the concept of transient energy growth. In §3, we present 
numerical results obtained using the eigenfunction approach for a variety of flows: plane 
Couette flow, constant angular momentum disk, Keplerian disk. In §4, we explain the 
physics of the numerical results by means of analytical and heuristic arguments. Finally, in 
§5, we discuss the implications of the results. In the Appendix, we describe the formalism 
to compute the transient energy growth. 

2. The Model 

2.1. Equilibrium Flow 

We consider a small patch of an accreting disk centered on radius Tq and viewed in a 
frame orbiting at the angular velocity of the gas at this radius. We employ Cartesian 
coordinates {X, Y, Z) such that X = r — ro is in the radial direction, Y — ro(0 — 0o) is in 
the azimuthal direction and Z is in the vertical direction. 

For ease of comparison with classical results in the fluid literature, we assume that the 
flow is incompressible, that it extends from X — —L to +L, and that there are rigid walls 
at the two ends with no-slip boundary conditions. The flow is unbounded along Y and Z. 
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In the limit L ro, the unperturbed velocity corresponds to a linear shear of the form 

V= (0,-^,0), (1) 

where Uo is the speed at the two walls. Because of rotation, a Coriolis acceleration acts on 
the fluid and is described by a frequency 

a;=(0,0,Qo), no^^. (2) 

Here the parameter q is positive (corresponding to angular velocity decreasing with 
increasing radius in a disk) and describes the radial dependence of fl{r) in the accretion 
disk, 

n{r) = Qo[jy. (3) 

Thus, q — 3/2 corresponds to a Keplerian disk and q — 2 corresponds to a disk with a 
constant specific angular momentum. For completeness, we note that q — 1 corresponds to 
a system with a flat rotation curve and g = to solid body rotation. 

The classical plane Couette flow that is widely discussed in the fluid hterature has a 
finite shear but no Coriolis force. In our model, this corresponds to a finite Uq but zero Qq, 
i.e., it represents the limit q — ^ cxd. The accretion disk problem, which is of primary interest 

in astrophysics, corresponds to finite q in the range 3/2 to 2. In comparing the present 
work to the fiuid literature, the reader is warned that our radial coordinate X maps to Y in 
the fiuid work, while our Y is their X. The notation wc use is standard in the astrophysics 
literature. Below we describe the self-contained set of generalized equations from beginning, 
for the convenience of the general reader. 



2.2. Perturbations 

The dynamics of a viscous incompressible fluid are described by the Navier-Stokes 
equation (e.g.. Landau & Lifshitz 1989), 

— + V.VV + Cu xuj X D + 2uj xV + V [-] = uV^V, (4) 
dt' \pj 

supplemented with the condition of incompressibility, 

V.V = 0, (5) 
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where t' is time, V is the velocity, Q is the Coriohs vector defined in equation (2), v is the 
kinematic coefficient of viscosity, D = {X,Y,Z), V' = {d/dX,d/dY ,d/dZ), and P is the 
pressure. Due to the incompressibihty assumption, the density p is a constant. 

It is convenient to analyse the perturbations in terms of dimensionless variables, 
X, y, t, defined by 

X = xL, Y = yL, Z = zL, V = UUo, t' = tL/Uo, (6) 

— * 

where C/ is a dimensionless velocity 

U={0,Uy,0), Uy^U{x)^-x. (7) 
Then, by substituting (6) into (4), we obtain 

— + U.WU + + + Vp = -V^f/, 8 

at q R 



where pU^ = P/p, d = {x,y,z), V = {d/dx,d/dy,d/dz), and the Reynolds number R is 
defined by 

ii=M. (9) 



We consider small perturbations in the velocity components of the form: 

Ux — u{x,y, z,t), Uy — > U{x) + v{x,y, z.t), Uz w{x,y, z,t), and perturbations in 
the pressure p ^ p + p{x, y, z, t). The linearized Navier-Stokes and continuity equation for 
the incompressible fluid then give 

fd 2v dp 

(^ + U—\v + u— + — + ^--V\ (U) 
\dt dy J dx q dy R ^ 

(I- -I) -1 4-^- (-) 

du dv dw ^ 
dx dy dz 

Let us rewrite the equations in terms of the x component of the vorticity, 

^ - ^ _ ^ (14) 
dy dz 
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Combining equations (10)-(13) and simplifying, we obtain 



^2 2 (dv du\ 

^ dx dy q \dx dy J 

Eliminating p and v from (10) by use of (13)-(15) we find 

Finally, combining (11) and (12) by use of (14) we obtain 

(^ + U— ](-——--—- -V^C (17) 
\dt dy J dx dz qdz R 

where we recall from equation (7) that U — —x. 

Equations (16) and (17) are the standard Orr-Sommerfeld and Squire equations, 
respectively, except that they now have additional terms proportional to 2/q because of the 
inclusion of Coriolis acceleration. We are interested in solving these linear equations with 
no-slip boundary conditions, i.e., u — v — w — QaX, the two walls. Equivalently 

du 

= — = C = 0, atx = ±l. (18) 
dx 

Because of translation-invariance of the unperturbed flow in y and 2;, we can decompose 
the perturbations in terms of Fourier modes in these directions. Also, for convenience, we 
study the perturbations in terms of (m, C) rather than {u,v,w). Therefore, we write the 
perturbations as 

— * 

u{x, y, 2;, t) — u{x, t) exp[iA;.r^] , 

C(x, y, z, t) = C{x, t) exp[ik.fp\, (19) 

where fp [= {y, z)] is any radius vector in the y — z plane and k = {ky, kz). By substituting 
(19) into (16) and (17), we obtain 

dC 

^ = -i[{jCc + Ccor)u + jC,gC], (20) 

where 

Cos = -(^' - kY'[{D^ - k'^f/iiR) - kyU{D'' - e) + kyD^U], 



-9- 



£c = -KDU, 
Csq = kyU-{D^-e)/{iR), 



^2 ,2\-l 



c — —C (n^ — k^ 

*^cor — ^cor\-'-^ 

D = d/dx. (21) 



If we further define 



Q=("), C=( ) , (22) 

\ C / \ ~l~ ^cor ^sq J 



equation (20) reduces to the form 



^ = -^CQ, (23) 



which we need to solve to obtain the eigenvectors and corresponding eigenvalues. Since the 
set of eigenmodes for this bounded flow problem is discrete and complete, we can write the 
solution to (23) in terms of an eigenfunction expansion, 

oo 

Q{x, t) = ^ ^Aj exp{-iXjt)Q^j{x) + Bj exp{-iij,jt)Q^j{x) 



u^j{x) 
u^j{x) 



where {Xj,Q^{x)) is the Orr-Sommerfeld eigensystem*^ and {iJ,j,Q'j{x)) is the Squire 
eigensystem. Formally merging the two systems, we can rewrite (24) as 

oo 

Q(x,t) = '^Cjexp{-iajt)Qj(x), 

- ( 11:)' ) . p^) 

where half the indices j correspond to the Orr-Sommerfeld modes and the other half of j 
correspond to the Squire modes, and aj — URj + iuij. Therefore, for the jth mode, (23) 
reduces to 

CQj = GjQj. (26) 



A complete set of eigenvalues and eigenvectors is called the eigensystem. 
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To calculate the set of eigenvalues and eigenvectors, we convert the differential operator 
L into an N X N matrix in a finite-difference representation and we then compute the 
eigenvalues and eigenvectors of the matrix. The required order N of the matrix for adequate 
accuracy depends on the physical parameters of the problem (mainly R and also ky, kz)- 
For the calculations presented here, we used N in the range 200 — 300 (i.e., each of the 
blocks Cos-i ^sq, and Ccor had a size in the range 100 to 150). 



The eigenvalues aj and the corresponding eigenvectors Qj for plane Couette flow have 
been studied by a number authors (e.g. Orszag 1971; Romanov 1973; Farrell 1988; Reddy 
& Henningson 1993) who have shown that there are no exponential growing eigenmodes in 
the system. That is, for no choice of the parameters is there an eigenvalue with positive 
(7/. Figure la shows such a typical eigenspectrum: the case shown has R — 2000 and 
ky — kz — 1- For comparison Fig. lb shows the eigenspectrum of a disk with constant 
angular momentum q — 2, and Fig. Ic a Keplerian disk q — 3/2. To the best of our 
knowledge, this one to one comparison of eigenspectra between standard plane Couette 
flow, a constant angular momentum flow and a Keplerian flow has not been reported earlier. 
It is clear that none of these flows has any growing eigenmodc and so all three systems are 
linearly stable. It should also be noticed that the eigenspectrum for plane Couette flow 
is very similar to that of a constant angular momentum flow. We explore this similarity 
further in §§3,4. 

One of the most important features of plane Couette flow (and its governing linear 
operator) is that the eigenmodes of the system are close to linearly dependent in nature i.e. 
they are non-normal in nature. Because of this, even in the absence of any exponentially 
growing mode in the system, it is possible to have a large transient growth in the energy of 
certain perturbations (Butler & Farrell 1992; Trefethen et al. 1993). This growth occurs in 
the absence of non-linear effects and is believed to play an important role in the transition 
from laminar to turbulent flow. 

Following previous authors (e.g. Trefethen et al. 1993; Schmid & Henningson 1994) we 
define the perturbed energy density as 



where a = 27: /ky, b = 271 /k^, and V = 2ab is the integration volume. We then seek to 
maximize the growth in this quantity. Recalling the formal solution of (23) in matrix form. 



2.3. Energy Growth 




(27) 



Q{x, t) — exp[—iCt\Q{x, 0), 



(28) 
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the maximum growth in the perturbed energy can be expressed as 



G{ky, kz, R, t) — maximum 



exp[-iCt]\\l, 



(29) 



where ||...||2 signifies the norm of the respective quantity and the subscript 2 specifies the 
2-norm or Euchdian norm. The 2-norm of the matrix can be evahiated by means of a 
singular vahie decomposition. Then, for a given t, the square of the highest singular value 
is the maximum energy growth, Gmax{t), for that time. Physically, by "maximum" we 
mean that we consider all possible initial perturbations Q{x, 0) and choose that function 
that maximizes the growth of energy at time t. The corresponding energy growth factor is 
G{ky, kz, R, t). When i = 0, by definition G{ky, k^, R, t) — 1, implying no growth. For given 
we maximize G{ky, k^, R, t) by writing Q{x, 0) as a linear combination of the eigenmodes 
of the system as in equation (25) and optimizing the coefficients Cj. The details are given 
in the Appendix. 

We should mention that to evaluate the growth one does not need to include all the 

eigenmodes in the computation. It has been shown by Reddy & Henningson (1993) that 
only a limited number, K/2, of the Orr-Sommerfeld and Squire modes, viz., those with the 
largest (i.e., least negative) aj values, are responsible for the growth. The remaining modes 
decay too rapidly to provide much growth. Therefore (25) can be rewritten as 



and the corresponding growth is GK{ky,kz,R,t). Here Q and C are N x K and K x 1 
matrices respectively and J^k is a K x K diagonal matrix consisting of the top K eigenvalues 
(top K/2 of Orr-Sommerfeld and K/2 of Squire eigenvalues). For the calculations presented 
in this paper, we generally used K < 60. 

The growth G{ky, kz, R,t) defined above is a function of four parameters. In various 
places in the paper we consider different kinds of maxima of this function. For instance, 
for fixed ky,kz, R, we could maximize G with respect to time t, and thereby determine 
the maximum growth G^axiky, kz,R)- This is the quantity that is plotted as contours in 
Figures 2 and 4. Or, we may wish to hold one of the components of the wavevector fixed, 
e.g., ky — (see Table 1, §4.1, also Figure 6 for other values of ky) or kz — (§4.2), and 
optimize the growth with respect to the other component of the wavevector and the time; 
this gives maximum growth factors such as Gmax{ky = 0, R) and Gmaxik'z = 0, /?). Finally, 
for a given i?, we could optimize over all the other parameters to calculate G„iax{R)- This 
is the quantity of most interest, and is shown for instance in Tables 1,2, and Figures 3, 5. 



K 



QK{x,t) = C-j exp{—iajt)Q-j{x) = Q exp[—iT,Kt]C 



(30) 
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3. Numerical Results 

3.1. Plane Couette Flow and Constant Specific Angular Momentum Flow 

In the previous section, we showed that the eigenspectra of plane Couette flow and 
a constant angular momentum disk (g = 2) are very similar. Here we show that the 
maximum growths are also similar. Figures 2a,b and 2c,d show contours of constant Gmax 
in the {ky,kz} plane for plane Couette flow and a g = 2 disk respectively for two values of 
R (500,2000). The maximum growth values for these two cases and for other values of R 
are found in Fig. 3. We see that the values are very similar for the two flows. Whereas 
for g = 2 the maximum growth occurs exactly on the kz axis [ky = 0), for plane Couette 
flow it is slightly off the axis, though by a progressively smaller amount with increasing 
R. This deviation in location of the occurrence of maximum growth in the ky — kz plane 
for a constant angular momentum disk compared to plane Couette flow was completely 
unnoticed earher, to the best of our knowledge. In cases of a large R, the best growth values 
for ky — compared to /Cj^ ~ (but ^ 0) do not have any practical difference. But for a 
small R the difference is important as in this case the maximum growth factor itself is small 
(see Fig. 3). Therefore, in presence of finite molecular viscosity this result has a physical 
implication to fiuid dynamics, though in the accretion disk when R is always expected to 
be very large this may not be an important issue. 

In Figure 3, we show the variation of the maximum growth Gmax and the corresponding 
time tmax at which the maximum growth occurs as functions of R for plane Couette flow 
and q = 2. We see that Gmax varies as R^ and tmax as R in both cases, with very similar 
values, again indicating the similarity of the two flows. However for plane Couette flow, 
growth maximizes for kz ~ 1.6, while for a g = 2 disk it happens at kz — 1.66. Moreover, in 
plane Couette Flow, the optimum ky scales as 1/R, already noticed by earlier authors (e.g. 
Butler & Farrell 1992), whereas for a g = 2 flow the optimum ky — 0. 

In the case of a constant angular momentum disk the epicyclic frequency of the disk 
becomes zero which makes the basic structure of the system very similar to that of plane 
Couette flow (compare the equation set (10)-(12) as well as (16) and (17) for plane Couette 
flow and constant angular momentum flow). This was already noticed by Balbus, Hawley & 
Stone (1996) and Hawley, Balbus & Winters (1999), who found from numerical simulations 
that these two flows are equally susceptible to hydrodynamic turbulence. We explore the 
physics of this similarity further in §4. 




k k 

y y 

Fig. 2. — Contours of Gmax{ky, kz, R) in the ky — kz plane, (a) Plane Couette flow for 
it! = 500: dashed contours correspond to G^ax — 2,4, ...,10, solid contours to Gmax — 
20, 30, 100, dotted contours to G^ax — 120, 140, 240, and dot-dashed contours to 
Gmax — 250, 260, 290. (b) Plane Couette flow for R — 2000: dashed contours correspond 
to Gmax = 10, 30, 130, solid contours to Gmax = 200. 300, .... 1000, dotted contours to 
Gmax = 1200, 1600, ...,4000, and dot-dashed contours to Gmax = 4500,4600,4700. (c) 
Constant angular momentum disk {q = 2) for R = 500: dashed contours correspond to 
Gm.ax = 2,4, ...,10, solid contours to Gmax = 20, 30, 100, dotted contours to Gmax = 
120, 140, 240, and dot-dashed contours to Gmax = 250, 260, 290. (d) Constant angular 
momentum disk for R — 2000: dashed contours correspond to Gmax — 10,30, 130, solid 
contours to Gmax = 200,300, ...,1000, dotted contours to Gmax = 1200, 1600, 4000, and 
dot-dashed contours to Gmax — 4500, 4600. 
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Fig. 3. — Gmax{R) and tmax{R) as functions of Reynolds number R. Solid lines correspond 
to plane Couette flow and dotted lines to a constant angular momentum disk {q = 2). The 
dashed hues show the analytic result discussed in §4.1.1. 
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3.2. Effect of Non-Zero Epicyclic Frequency 

The epicyclic frequency is given by 

K = ^2(2 - qn (31) 

which is zero for q = 2 and is non-zero (and real) for any q < 2. Figures 2(c),(d) show 
that, when q = 2, the maximum growth occurs for axisymmetric perturbations with vertical 
structure (hereafter we refer as "vertical perturbation"), i.e. ky = 0, fc^ 7^ 0. We begin this 
section by exploring what happens to such vertical perturbations when q < 2. Note that 
most of the fluid literature is devoted to plane Couette flow with no Coriohs force. In the 
astrophysics literature, Yecko (2004) mainly concentrated on the Keplerian disk (q — 1.5) 
and plane Couette flow, whereas we consider the full range, 1.5 < q < 2. 

Table 1 shows the maximum growth for vertical perturbations for four values of q: 
1.99, 1.9, 1.7, 1.5. We see that, as q decreases, the growth falls dramatically, and so does 
the time at which the maximum occurs. For a Keplerian flow, the growth factor is under 
4. Moreover, numerical experiments show that the growth is insensitive to the value of R. 
In other words, the growth is not limited by viscosity but rather by the dynamics itself. 
Wc explore the reasons for this in §4.1. Similar results are discussed in AMN05, using a 
Lagrangian picture. 

Table 1 

Energy Growth Factors of Disks for R = 2000 and ky = 



Q 


kz 


Gmax{ky^O,R^ 2000) 


tmax{ky^O,R^ 2000) 


1.5 


2.5 


3.84 


2.8 


1.7 


2.1 


6.31 


1.1 


1.9 


2.1 


18.16 


8.3 


1.99 


1.9 


153.4 


27.8 


2 


1.66 


4661 


277 



We next remove the restriction to vertical perturbations and consider general ky, k^- 
Figure 4a shows contours of constant growth for q = 1.99 for R = 2000. Even though the 
value of q is only very slightly different form that used in Figure 2d (g = 2), nevertheless 
we sec a dramatic change. The main qualitative difference between the two cases is that 
the epicyclic frequency is zero for q — 2 but is (slightly) non-zero for q — 1.99. Still, 
this small changes causes a major modification of the results, showing what a dominant 
effect the epicychc frequency has on the fluid dynamics. The other panels in Fig. 4 show 
results for other values of g. It is interesting to see how the location of the maximum in 
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the ky — kz plane changes as the system approaches the Keplerian regime, and also how 
the magnitude of the growth reduces. This change in location of maximum growth as a 
function of q in the ky — k^ plane has not been noted earlier, to the best of our knowledge. 
For q = 2 the maximum energy growth is a factor of 4600 and occurs on the ky = line 
(Fig. 2d) while for g = 1.5 the maximum growth is only 22 and occurs on the kz = line 
(Fig. 4d). Therefore for a constant angular momentum disk, we need to include vertical 
structure in the perturbations to maximize the energy growth, whereas for a Keplerian disk 
a 2-dimensional analysis is sufficient. Table 2 lists the maximum energy growth factors for 
g = 1.5 - 2 when R = 2000. 

Figures 5a and 5b show respectively the variation of the maximum growth and the 
time at which the maximum growth occurs as a function of Reynolds number in a Keplerian 
disk. It is seen that for large i?, Gmax scales as R"^^^ and tmax as R^^^. This suggests that, 
even though the growth is modest for the values of R we have considered, if we go to 
sufficiently large values of R, very large energy growth might still be possible. This is of 
interest because the Reynolds number of a cold accretion disk is very high (many orders of 
magnitude higher than the values considered in this work), so turbulence could be generated 
in such systems. Due to numerical constraints our current results are limited to i? < 10"^; 
however, this range captures most of the basic features of the growth. Yecko (2004) used a 
superior spectral code and was able to go to much larger values of R. 

Figure 6 shows how Gmax and tmax scale with ky at a given R. The maximum growth is 
achieved at ky ~ 1.2. At smaller ky, Gmax scales as ky^^, while at larger ky, Gmax decreases 
as ~ ^/ky or ky^^^. Also tmax scales as k'"^^^ at large t. Yecko (2004) and Umurhan & 
Regev (2004) identified the scaling of Gmax with R for a Keplerian disk, but the other 
scahng relations have not been discussed before. We also derive all the scahng relations 
analytically, for the first time, in §4. A detailed understanding of these scalings is given in 
§4.2. Identical scahng relations are also derived by AMN05. 
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Table 2 

Maximum Energy Growth for Disks with Various Values of q and R — 2000 



<1 


ky 


kz 


GmaxiR = 2000) 


tmax{R = 2000) 


1.5 


1.2 





21.67 


11 


1.7 


1.05 


0.6 


23 


12.4 


1.9 


0.34 


0.96 


32.33 


19.8 


1.99 


0.11 


1.61 


171.13 


31.3 


2 


C) 


1.66 


1661 


277 



3.3. Nature of the Growing Perturbations in a Keplerian Flow 

We have seen that for a Keplerian flow the maximum growth occurs for ky ~ 1.2, 
kz — 0. Figures 7 and 8 show the development with time of the perturbed velocity 
component u{x,y) corresponding to R = 500,4000, respectively, optimized for the best 
growing mode. In each case, we show snapshots corresponding to four times: t = 0, tmaxf^, 
tmax, ^tmaxf^- The perturbatious are seen to resemble plane waves that are frozen in the 
shearing flow. The initial perturbation at f = is a leading wave with negative kx and 
with \kx\ ^ ky. With time, the wavefronts are straightened out by the shear, until at 
t — tmax, the wavefronts are almost radial and k^ ^ 0. At yet later times, the wave becomes 
trailing and the energy also decreases. The perturbations are very similar to the growing 
perturbation described by Chagelishvili et al. (2003) and Umurhan & Regev (2004). 
However, those authors (and also AMN05) considered an infinite system whereas our fluid 
is conflned to a box of size 2L in the x direction. Figure 9 shows the optimum growth of 
the energy G{t) as a function of time for the two perturbations whose time evolution for 
the best growing modes are shown in Figures 7 and 8. 



4. Physical Interpretation of the Numerical Results 

In this section we attempt to understand via an analytical approach the numerical 
results of the previous sections. We also like to derive the scaling relations described in the 
previous sections analytically. We show that the analytical solutions match the numerical 
results well. In the interest of clarity, we work with the original dimensioned equations. 
Thus X (going from —L to +L), Y, Z are our coordinates, and we write the corresponding 
components of the wavevector as kx, ky, kz, respectively. We use t for time (called f in 
§2.1). 
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Fig. 4. — Contours of Gmaxi^y, k^, R) for disks with different values of q and R = 2000. 
(a) q — 1.99: dashed contours correspond to Gmax = 2,4, ...,10, sohd contours to 
Gmax — 15, 30, 105, dotted contours to Gmax — 130, 140, 170, and dot-dashed contours 
to Gmax = 174, 174.2, 174.4. (b) q = 1.9: dashed contours correspond to Gmax = 2,4, 14, 
sohd contours to Gmax — 16, 19, 28, and dotted contours to Gmax — 29, 29.5, 32. 
(c) q = 1.7: dotted contours correspond to Gmax — 2,3, ...,9, and sohd contours to 
Gmax = 11, 12.2, 23. (d) q = 1.5: dashed contours correspond to G,nax = 2,3, ...,6, 
sohd contours to Gmax — 7, 9, 13, and dotted contours to Gmax — 17.5, 18, 21.5. 
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Fig. 5. — Solid lines show (a) Gmax{R), and (b) tmax{R)-, as functions of it! for g = 1.5. 
Dotted lines correspond to the analytical result for kx,min — 1-7 discussed in §4.2. 



- 21 - 




-22- 




c 



3 



2 









X 




Fig. 7. — Shows the development of the perturbed velocity u{x^ y) as a function of time 
for the best growing energy optimal in a Keplerian flow with R — 500. The perturbation 
has ky — 1.29, — 0, and the maximum growth is achieved at tmax{R — 500) = 6.6. 
The four panels correspond to (a) t = 0, (b) i = tmax/'^ = 3.3, (c) t = tmax = 6.6, (d) 
t = Stmaxf^ = 9.9. Solid and dotted contours correspond to positive and negative values of 
u respectively. 
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Fig. 8.— Same as Fig. 7 but for R = 4000. Here ky = 1.2, = 0, tmax{R = 4000) = 13.3, 
and the four panels correspond to (a) t — 0, (b) t — tmax/'^ — 6.65, (c) t — tmax — 13.3, (d) 
t — St max /"^ = 19.95. Solid and dotted contours correspond to positive and negative values 
of u respectively. 
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Fig. 9. — Growth of the perturbed energy G{t) as a function of time for the cases shown in 
Figs. 7 and 8. 



-25- 



Let us define the shear frequency 2 A and the vorticity frequency 2S as follows 
(Narayan, Goldreich & Goodman 1987), 

2A^-qn, 2B = 2{A + n) = {2-q)n. (32) 

Then the three components of the momentum equation and the incompressibility condition 
give 



du ^ dp „,2 


(33) 


dv „ dp „,2 


(34) 


dw dp ,2 


(35) 


du dv dw 
dX^W^dZ^ ' 


(36) 



where p = P/p. The Lagrangian time derivative d / d,t is given by 

The numerical results described in §3 showed that plane Couette flow and q — 2 flow 
both have maximum energy growth for perturbations with vertical structure, whereas 
Keplerian g = 3/2 flow has maximum growth for two-dimensional perturbations with no 
vertical structure. We analyse these two cases separately 



4.1. Vertical Perturbations 

To begin with, let us ignore the walls and assume plane wave solutions of the form 

u{X, F, Z, t) = u{t) exp{ikxX + ikyY + ikzZ), etc.. (38) 

Furthermore, since perturbations with ky — (equivalent to ky — 0) were seen to 
grow robustly for both plane Couette flow and q — 2, let us assume ky — 0. For such 
perturbations, d/dt — d/dt. 



4.1.1. Plane Couette Flow 



For plane Couette low, we set 2Q = in equation (33) and 2B — 2 A in equation (34). 
We can show that the fastest growing plane wave perturbation for a given kx, kz takes the 



-26- 



form 



u — uq exp[ikxX + ikzZ — ^{k^ + 
V — —2AtuQ ex:p[ikxX + ikzZ — ^{k^ + 

kx 

-uq exp[ikxX + ikzZ - v{k\ + k%)t\, 



w 



k' 



(39) 
(40) 

(41) 

(42) 



where uo is an arbitrary amplitude. The ratio of energy at time t to the initial energy is 
then given by 



G{t) 



7/2(0) +^2(0) +^2(0) 



1 + 



A-2 

^ {2AtY 



exp 



-2u{kx + kl)t 



(43) 



Since we are interested in flows that have perturbations with large growth, let us ignore the 
1 in the first factor. Then, the time at which the energy is maximum is given by 



t. 



y{k\ + kiy 
and the corresponding energy growth factor is 



(44) 



2A\ 



(45) 



The problem we have analysed in §3 is a flow with walls at X — ±L with no-slip 
boundary conditions. In the absence of viscosity, the simplest solution to this problem is 



kz ( 'kX\ 

u — —- — w — Uq (1 + cos —— exjpnkzZ), 
kx V L ' 



V — p — 0, 



(46) 



which can be seen by inspection to satisfy the boundary conditions at X = ±L. This 
solution is the sum of a plane wave with kx = and amplitude Uq, and two waves with 
kx — ±7r/L and amplitude uo/2. Roughly, we expect that the solutions for tmax, Gmax we 
wrote down earlier for a single plane wave will be approximately correct provided we set 
kx equal to (7r/L)2/2, i.e., the mean of and {n/Ly. Noting that the Reynolds number is 
given by 

(47) 



V 

we find that the maximum growth factor is given by 

Gmax{kz, R) 



R^klL^ e-2 



|7r2 + A;2 L2 



(48) 
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Maximizing this over kz-, we find that the optimum wavevector is 

kzL = k, = Ti/2 = 1.57. (49) 

This is close to the numerically determined value of 1.6 given in Table 1. The maximum 
growth factor and time at which this maximum is attained are then 

GmaxiR) = 0.82 X lO-^i?^ (50) 

\2A\tmax{R) = 0.13i?. (51) 

These relations are shown as dashed fines in Figure 3. We see that the scaling with R agrees 
well with the numerical results, and the coefficient is also reasonably close. 



4- 1-2. Constant Specific Angular Momentum Flow 



In this section we consider a rotating flow with q = 2. The vorticity frequency 2B 
vanishes, and so the term proportional to it is not present in equation (34). We can then 
write down the following plane wave solution 



u — 



k\ + k% 



2ntvo exp[ikxX + ikzZ - iy{kx + kl)t], 



w 



V — Vq exp[ikxX + ikzZ — i'{k\ + k\)t\^ 
kxkz 



kx + kz 

ikx 



where vq is an arbitrary amplitude. 



2fltvo exp[ikxX + ikzZ — u^k^ + k\)t], 
2VLvq ex.\)[ikxX + ikzZ — v{k\ + A;|)t], 



(52) 

(53) 

(54) 

(55) 



This solution looks quite different from that for plane Couette fiow. For instance, here 
u and w grow linearly with time at early time and v remains constant, which is the reverse 
of the case for plane Couette flow. Also, now we have a non-zero pressure perturbation. 
Nevertheless, the energy growth factor has the same dependence on time as for plane 
Couette flow: 



G{t) 



m2(0) + 't;2(0) + w2(0) 



1 + 



kx + kz 



{2Vttf 



exp 



-2v{k\ + kl)t 



(56) 



Note that, for q — 2, 20, — 2A, so the result is in fact identical. The reason for this 
close similarity is apparent when one considers the original dynamical equations (33)-(35). 
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The only difference between plane Couette flow and q — 2 flow is that for the former the 
term 2^^ is missing in equation (33) and the term —2Bu present and is equal to —2Au 
in equation (34), whereas for the latter the term 2flv is present and is equal to 2Av in 
equation (33) and the term —2Bu is missing in equation (34). The equations are thus very 
symmetrical, except that X and Y are interchanged in the two cases. The resuhing flows 
look very difi^crent because of the switch in coordinates, but the growth is identical 

The rest of the analysis proceeds exactly as in the previous subsection. As before, 
we conclude that the optimum kzL = ~ 1.57, that the maximum growth factor is 
Gmax ~ 0.82 X 10~^i?^, and that the maximum growth happens at a time tmax ~ 0.13i?. As 
we saw in §3, the numerical results are indeed very similar for plane Couette flow and q = 2 
flow (Figure 3). The present analysis explains why this happens even though the dynamics 
are quite different. 



Now we consider a more general flow with q < 2. Such a flow has an angular 
momentum gradient that is stable according to the Rayleigh criterion. This leads to 
epicyclic oscillations with frequency k [see Eq. (31)]. 

Let us ignore viscosity. Then, deflning k = {kx + /^|)^^^, the plane wave solution with 
maximum growth is given by 



4.1.3. q<2 Flow 




(57) 




(58) 



(59) 




(60) 



The energy growth as a function of time is given by 




(61) 



Clearly, the maximum possible growth is 



G. 



2 



(62) 



max 




-29- 



and the growth happens after a time proportional to one quarter the epicychc period, 



In the actual flow with walls and viscosity the growth will be a little less, as confirmed 
by the numerical results in Table 1, but equation (62) gives a rigorous upper limit to the 
growth. 

The key point of this analysis is that, for q < 2, there is a limit to the growth that arises 
just from dynamics; specifically, it is caused by the presence of a non-vanishing epicyclic 
frequency. Moreover, the limiting growth in energy is just a factor of 4 for a Keplerian fiow. 
The hmit has nothing to do with viscosity. In contrast, both plane Couette fiow and q — 2 
fiow can have infinite growth as far as the dynamics are concerned and the limit to growth 
arises only from viscosity. 

Balbus, Hawley & Stone (1996) and Hawley, Balbus & Winters (1999) suggested that 
the existence of epicyclic motion lends dynamical stability to fiows with q < 2 and that 
this makes these flows more resistant to turbulence. Our analysis confirms their suggestion 
for perturbations with vertical structure. However, their argument does not apply to the 
two-dimensional perturbations we consider next. 



The last subsection explained why perturbations with vertical structure have very little 
growth for q < 2. The growth is especially insignificant for a Keplerian flow. §3 showed that 
these fiows have more growth for perturbations with non-zero ky- In fact, for a Keplerian 
fiow, the maximum growth is for kz = 0, ky ^ 0. We now consider such perturbations. 

We consider a plane wave that is frozen into the fiuid and is sheared along with the 
background flow (see Figs. 7, 8). If the fiow starts at time t — with initial wave- vector 
{kxi, ky) in the XY-plane, then the X-wavevector at later times is given by 




(63) 



4.2. 



Two-Dimensional Perturbations 



kxit) — kxi + q^kyt. 



(64) 



With the above definition of kx, we consider a plane wave solution of the form 



u — u{t) exp{ikxX + ikyY), etc.. 



(65) 



Because of the non-zero ky, the Lagrangian time derivative is given by 



(66) 
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The relevant plane wave solution in the absence of viscosity has been written down by 
a number of authors (e.g., Chagelishvili et al. 2003; Umurhan & Regev 2004). Generalizing 
the solution for finite viscosity, we have 



k 

u 



C-^ exp (^ikxX + ikyY -V k\t')dt'^ , (67) 

V = -C^ exp (^ikxX + ikyY -u k\t')dt'^ , (68) 
w^O, (69) 
p^i( (^-^2n - ^2gflj exp (^ikxX + ikyY - v j'^ k\t')dt'^ , (70) 

where C is the amplitude of the vorticity perturbation. Since w = 0, we see that the 
perturbations are two-dimensional (hence the name "two-dimensional perturbations"). 
Also, the velocity components are independent of q and it is the pressure that adjusts so as 
to keep the dynamics the same for all values of q. In fact, the above solution is valid even 
for plane Couette flow, provided we make the replacements 20 ^ and 2qfl —2A. 

In the absence of viscosity, the energy growth is given by 

rox ' n,y 

that is, the energy is inversely proportional to the square of the total wave-vector 
fc^ = k\ + ky. This result is easy to understand. For inviscid incompressible two- 
dimensional flow, the vorticity V x is exactly conserved. This means that k-d is constant, 
so the velocity scales inversely as k. The energy must then vary as A;~^. The energy is thus 
largest when k is smallest. Using equation (64) we now see what is required if we wish to 
obtain a large energy growth. We need to start with a large negative value for kxi- As time 
goes on, kx will become progressively less negative; as a result, k will decrease and G will 
increase. The maximum growth will be achieved when kx — ^i giving 



^ \0A\f 

Ky Ky 



Now consider the effect of having rigid walls at X = ±L. By the uncertainty principle, 
kx cannot become exactly zero, but must have a minimum magnitude, kx,min ~ tt/L. The 
maximum energy growth is then approximately given by 



^ = = kCin+'k'y ^ (W + W ^^^^ 
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where we have assumed that kxi ^ ky- Including also the effect of viscosity this becomes 

.„«)._il|_exp(-A^I^). p4) 

Maximizing this with respect to kxi and /cy, we obtain 

kyL^ky^ 2.2, kxiL = k^i ~ l.m^'\ (75) 

V2 

The maximum growth and the time of maximum are then 

Gmax{R) ~ 0.059i?2/3^ \2A\tmax{R) ~ 0.59/?^/^ (76) 

While the scalings with R are accurate and agree with the numerical results presented 
in §3, the coefficients are approximate since they depend on the assumed value of 
kx,min- If instead of taking kx,minL = vr, we select kx,minL = 1, we find kyL = 0.71, 
kxiL = 0.89i^l/^ G^axiR) ~ 0.27i?2/3^ \2A\t^ax{R) ~ l.25R^/^. The results in §3 lie 
between these two estimates. In fact, if we choose kx. minL ~ 1-7, then we obtain kyL — 1-2, 
Gmax{R) ~ 0.13i?2/3, \2A\tmax{R) ~ O.SSi?^/^, which agree with the numerical result (Table 
2). 

Finally, we can carry out the analysis by keeping ky fixed and optimizing only kxi- We 
then find that Gmax{ky, R) and tmax{ky, R) vary as 

\2A\t^,,{ky, R) = {kyL)-^/^R^/\ (78) 

We see that Gmax varies as (/cyL)^/^ for small kyL and as {kyL)~^^^ for large kyL. The 
time of maximum scales as [kyL)'"^^^ for all kyL. The above analytical results are plotted 
as dotted lines in Figure 6, assuming kx,minL — 1.7 as derived above. We see that the 
agreement with the numerical results is very good except at very large R where the 
calculations are not very accurate. 



5. Discussion and Conclusions 

We have demonstrated that significant transient growth of perturbations is possible 
in a Keplerian flow between walls (as shown by Yecko 2004). Although the system does 
not have any unstable eigenmodes, nevertheless, because of the non-normal nature of 
the eigenmodes a signiflcant level of transient energy growth is possible for appropriate 
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choice of initial conditions. If the maximum growth exceeds the threshold for inducing 
turbulence, it is plausible that this mechanism could drive the system to a turbulent state. 
Presumably, once the system becomes turbulent it can remain turbulent as a result of 
nonlinear interactions and feedback among the perturbations. 

In this so-called bypass mechanism for transition to turbulence, the maximum energy 
growth and the time needed for this growth are likely to be the main factors that control 
the transition to hydrodynamic turbulence. It has been observed in laboratory experiments 
that plane Couette flow can be made turbulent for Reynolds numbers above a critical 
value Rc ~ 350. According to our analysis, for R = 350, the maximum energy growth is 
Gmax{R = 350) = 145, and the maximum occurs at time tmax = 42.3 (Fig. 3). Since a 
constant angular momentum disk {q = 2) behaves very similarly to plane Couette flow, the 
critical Reynolds number for turbulence for this case is also likely to be Rc ~ 350. For this 
R, the growth factor is Gmax{R — 350) = 143.5 and the time-scale is tmax — 48.3. 

It is true that the underlying equations as well as the maximally growing modes in 
plane Couette flow (and a q—2 disk) are different compared to a Keplerian flow, and so 
the turbulent phases in the two systems may have significant differences. However, the 
presence of similar boundary conditions may suggest similarities in the kinematic structure 
of the turbulence. In fact, a Keplerian disk and a constant angular momentum disk are 
two special cases of rotating shearing flows parameterized by q. Based on this, we make 
the plausible assumption that the threshold energy growth factor needed for transition to 
turbulence in a shear flow with any value of g is i?c ~ 145. Applying this conjecture to the 
optimal two-dimensional perturbations of a Keplerian disk analysed in §4.2, we estimate 
the critical Reynolds number for a Keplerian flow to be Rc ~ 3.4 x 10"^, i.e., a factor of 100 
greater than for plane Couette flow. The time to reach the maximum is t^ax — 28.3, which 
is comparable to that in plane Couette flow, and is not too large compared to the accretion 
time scale of a geometrically thin disk. 

Instead of taking Rc ~ 350, which is perhaps somewhat optimistic since plane Couette 
flow needs to be perturbed signiflcantly before it will become turbulent at this Reynolds 
number, we might wish to be conservative and assume Rc ~ 1000 for this flow. At this value 
of R, plane Couette flow and q = 2 flow have Gmax{R = 1000) ~ 1200 and tmax ^ 120 — 140. 
Applying the requirement Ec ~ 1200 to Keplerian flow, we find Rc ~ 10^ and tmax ~ 100. 
Now the critical Reynolds number is a factor of 1000 greater than for plane Couette flow. 

Why is the critical Reynolds number so much larger for a Keplerian disk compared to 

a constant angular momentum disk or plane Couette flow? The numerical results in §3 
and the analytical work in §4 provide the answer, viz., the presence of epicyclic motions 
in a Keplerian disk. It is very interesting to note that the presence of epicyclic motion 
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not only kills growth dramatically, it also changes the optimum wavevector {ky,kz} of the 
perturbations needed to produce energy growth. For a constant angular momcntTim disk 
{q — 2) and plane Couette flow, both of which have zero epicyclic frequency, it is seen that 
growth is maximized for ky ^ (on the kz axis). Even for a very small shift in the value 
of q below 2, corresponding to the introduction of a small epicyclic frequency, the location 
of maximum growth immediately moves signiflcantly in the ky — k^ plane from the kz axis 
(see Fig. 4 which corresponds to R — 2000). With decreasing g, the epicychc motion of the 
disk increases, and correspondingly the optimum value of ky for growth increases while the 
optimum kz decreases. When q — 1.5, i.e., when the disk is exactly Keplerian, the growth 
is maximum for kz — (on the ky axis). To the best of our knowledge, this change in the 
location of the maximum in the ky — kz plane has not been commented upon prior to this 
work. 

The change between q — 2 and g = 1.5 may be completely understood analytically, 
as we show in §4. The important point is that the vertical perturbations {ky — 0) that 
cause the large observed growth in a g = 2 disk require an absence of epicyclic motions. 
When the epicyclic frequency is zero, the velocity perturbation is able to grow linearly with 
time and the energy grows quadratically. The only limit to growth is provided by viscosity, 
which gives a scaling Gmax oc R^. However, once there is a non-zero epicyclic frequency, 
the growth is immediately limited. Even in the absence of viscosity, only a modest level of 
growth is possible. In fact, for a Keplerian flow, the maximum growth that one can obtain 
from vertical perturbations is only 4, well below the critical growth needed for turbulence. 
If vertical perturbations were the sole route to turbulence, then a Keplerian flow could 
never make the transition to turbulence. 

However, as §4.2 shows, there are other kinds of perturbations, speciflcally two- 
dimensional perturbations with kz — 0, which are not affected by epicyclic motions. For 
these perturbations, pressure fluctuations are able to absorb the effect of the Coriolis force. 
As a result, two-dimensional perturbations are able to grow to arbitrarily large values in 
the absence of viscosity. However, the growth is much reduced compared to the vertical 
perturbations described in the previous paragraph and it scales only as i?^/"^. Thus, one 
needs much larger values of i? ~ 10^'^ — 10® before one can achieve the same level of energy 
growth as can be found in a g = 2 disk for Reynolds numbers as small as 10^'^ — 10^. 

These results lead to a better understanding of the numerical simulations described 
in Balbus, Hawley & Stone (1996) and Hawley, Balbus & Winters (1999). Both papers 
showed that there is a close similarity between plane Couette flow and q = 2 flow, in the 
sense that the two flows readily became turbulent in numerical simulations. However, once 
the authors reduced the value of q below about 1.95, no turbulence was seen even when the 
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flows were initialized with large perturbations. The authors suggested that the change in 
behavior is because of the dynamical stability imposed by the Coriolis force and epicyclic 
motions. Our analysis supports this conclusion. 

However, Balbus, Hawley & Stone (1996) and Hawley, Balbus & Winters (1999) then 
proceeded to rule out the possibility of hydrodynamic turbulence in Keplerian disks. We do 
not agree with this conclusion. As we have shown, Keplerian disks can indeed support large 
transient energy growth, but they need much larger Reynolds numbers to achieve the same 
energy growth as plane Couette flow or g = 2 flow. The numerical simulations probably 
had effective Reynolds numbers <^ 10^ (because of numerical viscosity) which is below our 
most optimistic estimate of the critical Reynolds number. Thus, we suspect the simulations 
did not have sufficient numerical resolution to permit turbulence. In fact, Longaretti (2002) 
already suspected that the non occurrence of turbulence in previous simulations may be 
just due to the choice of low Reynolds number. 

Although the problem we analysed is shear flow between walls, the optimum growing 
perturbations that we find for the Keplerian case are very similar to those described by 
Chagelishvih et al. (2003) and Umurhan & Regev (2004) for an infinite shear flow. The 
perturbations are basically plane waves that are frozen in the shearing flow. Initially, at 
t — 0, the effective wave vector of the perturbation in the x direction {kx) is negative, which 
means that we have very asymmetric leading waves. As time goes on, the wavefronts are 
straightened out by the shear and \kx\ decreases. At the time when the growth is maximum, 
kx ^ (but not precisely because of the walls, see §4.2) and the wavefronts become almost 
radial. At yet later time, the growth decreases and the wave becomes of a trailing pattern. 

The above time evolution is very different from that seen for the optimum perturbations 
in plane Couette ffow or in a g = 2 disk. In plane Couette ffow, the x-component of the 
perturbation, u (i.e. the normal velocity), dominates over the other components, v, w, at 
t = 0. However, u remains at the same level for all time whereas v and w increase strongly 
up to the point of maximum growth before declining. The overall shape of the perturbation 
is roughly self-similar with time. For a constant angular momentum disk, on the other 
hand, it is v which remains constant with time whereas u and w vary by large amounts. 
However, as in plane Couette flow, the solution is largely self-similar in character up to the 
maximum. Neither of these flows shows the shearing perturbations that are characteristic 
of the Keplerian problem (Figs. 7 and 8). 

We conclude with an important caveat. While the demonstration of large energy 
growth is an important step, it does not prove that Keplerian disks will necessarily become 
hydrodynamically turbulent. Umurhan & Regev (2004) have shown via two-dimensional 
simulations that chaotic motions can persist for a time much longer than the time scale 
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tmax needed for linear growth. However, they also note that their perturbations must 
ultimately decline to zero in the presence of viscosity. To overcome this limitation, it is 
necessary to invoke three-dimensional effects. Secondary instabilities of various kinds, 
such as the elliptical instability, are widely discussed as a possible route to self-sustained 
turbulence in linearly perturbed shear flows (see the review article by Kerswell (2002); see 
also e.g. Hellberg & Orszag 1988; Le Dizes & Rossi 1996). It remains to be seen if these 
three-dimensional instabilities are present in perturbed flows such as those shown in Figures 
7 and 8. If they are, one will in addition have to show that they lead to non-linear feedback 
and 3-dimensional turbulence. 

We would like to thank the referee for various suggestions that improved the 
presentation of the paper. This work was supported in part by NASA grant NAG 5- 10780 
and NSF grant AST 0307433. 



A. Appendix: Method to Compute the Transient Growth 

To compute the optimum growth, flrst we need to evaluate the 2-norm of Q. Prom (29) 
it is clear that the 2-norm depends on C which consists of Cos ^-nd Csq- The underlying 
Hilbert space of the second order linear operator, jCgq: is Hgq — -^^[—1, 1] ^- Therefore the 
inner product of Ci; C2 G is deflned as 

(Ci,C2)l = ^aci^^^- (Ai) 

The domain of Cgq, that is Vgg, is the set of functions {■0} which have a second derivative 
in 1, 1] satisfying ip{±l) = 0. Following DiPrima & Habetler (1969) we also deflne the 
underlying Hilbert space of Cos, 'Hos, consisting of the set of functions having a second 
derivative in 1, 1] satisfying '^(±1) = 0. Therefore, for Ui,U2 G Ti-os the inner product 
is deflned as 

{ui, U2)h = {Dui, Du2)l + k^{ui, U2)l- (A2) 

The domain of Cos, that is Vos, is the set of functions that have a fourth derivative in 
L^[— 1, 1] satisfying — ^''(il) — 0- Therefore the underlying Hilbert space of C is 

Ti = Tios X T^sq and the corresponding domain is I? = x P^y. Thus combining (Al) and 
(A2) and with some algebra, the inner product for Qi,Q2 & Ti- can be written as 

{Qi, Q2) = (^tii, U2)l + (Ci, C2)l, (A3) 



^The Hilbert space is defined as a complete vector space with an inner product. 
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where = -(^2 _ 

Now following Butler & Farrell (1992) the perturbation energy density can be evaluated 



as 



= -L /' r i\u^ + ^2 ^ w^)dzdydx, (A4) 
2V J -I Jo Jo 

where a = 2TT/ky, b = 2iT/kz, and V = 2ab is the integration volume. The physical velocity 
components are the real quantity obtained as 

1 . 

u = -{u exp[ik.fp] + u* exp[—ik.fp]}. (A5) 

Now replacing {v^ + w^) in terms of C and du/dx in (A4) and integrating over y and z we 
obtain 



F- ^ 



"''''^ dx dx^^^ 



dx, (A6) 



where u and ( are considered to be N dimensional column matrices. Now combining (30), 
(A3) and (A6) we obtain 

8k^E = WQkW^ = C^e^^^*Qe-^^^*C, (A7) 

where Q is a x Hermitian matrix whose ij'th element is the inner product of Qi and 
Qj, 

Qij = {Qi: Qj) = {^Uh Uj)L + {Ch Ql- (A8) 

Decomposing Q in terms of a matrix W according to Q = W'^W, and combining (29) and 
(A7) , we obtain the expression for the optimum growth (see also Reddy & Henningson 
1993 and Schmid & Henningson 1994), 

GKit) ^ maximum f l'^^^^P^'^.f]^''^ ^ = \\W exp[-iS;,t] W-'\\l (A9) 



\WC 



12 



where the subscript 2 denotes the 2-norm or Euclidian norm. The 2-norm of the matrix 
W exp[— iExi] can be evaluated by means of a singular value decomposition. Then, for 
a given t, the square of the highest singular value is the maximum energy growth, Gmax{t), 
for that time. W can be computed easily by a similarity transformation of Q 

Q = S\fQ'dS^S\/^dS^ = W^W, (AlO) 
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where S is a unitary matrix and Qa is a diagonal matrix consisting of the eigenvalues of Q 
along the diagonal. Therefore one only needs to construct the matrix Q to compute Gxit), 
while is immediately available from the eigenvalues of £. From (A6), (A7) and (A8), in 
the finite-difference approximation Q can be written as 



Q 



Ax 



ox ox 



■'mj 



mj 



(All) 



where m runs between 1 and N (points on the finite-difference grid) and j ranges from 1 to 
K (mode numbers). 
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